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ABSTRACT 

Context. The hydrogen ionisation degree deviates substantially from statistical equilibrium under the conditions of the solar chromosphere. A 
realistic description of this atmospheric layer thus must account for time-dependent non-equilibrium effects. 

Aims. Advancing the realism of numerical simulations of the solar chromosphere by improved numerical treatment of the relevant physics will 
provide more realistic models that are essential for interpretation of existing and future observations. 

Methods. An approximate method for solving the rate equations for the hydrogen populations was extended and implemented in the three- 
' — ^ dimensional radiation (magneto-)hydrodynamics code C05BOLD. The method is based on a model atom with six energy levels and fixed 
radiative rates. It has been tested extensively in one-dimensional simulations. The extended method has been used to create a three-dimensional 
model that extends from the upper convection zone to the chromosphere. 
O ' Results. The ionisation degree of hydrogen in our time-dependent simulation is comparable to the corresponding equilibrium value up to 
jf-H | 500 km above optical depth unity. Above this height, the non- equilibrium ionisation degree is fairly constant over time and space, and tends 
ry} , to be at a value set by hot propagating shock waves. The hydrogen level populations and electron density are much more constant than the 
corresponding values for statistical equilibrium, too. In contrast, the equilibrium ionisation degree varies by more than 20 orders of magnitude 
^ between hot, shocked regions and cool, non- shocked regions. 
• t-h Conclusions. The simulation shows for the first time in 3D that the chromospheric hydrogen ionisation degree and electron density cannot be 
1 calculated in equilibrium. Our simulation can provide realistic values of those quantities for detailed radiative transfer computations. 

u 
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1. Introduction 

Hydrogen is the most abundant constituent of the solar gas 
and the major electron donor in the solar chromosphere. The 
assumption of local thermal equilibrium (LTE) for hydrogen 
atoms is not valid in the chromosphere because the radiative 
transitions rates dominate over the collisional rates. Even the 
assumption of instantaneous statistical equilibrium, where the 
level populations of each species at all times are in equilib- 
rium as determined by the local thermodynamic state and the 
non-local radiation field, is not valid in the dynamical chromo- 
sphere. The timescale on which the hydrogen level populations 
adjust to changes in the atmosphere (which is set by the transi- 
tion rates between levels) is too long compared to the timescale 
on which the atmosphere changes. In order to properly model 
this behaviour one has to solve the rate equations for the popu- 
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lations of all relevant energy levels. We henceforth refer to this 
approach as time-dependent non-LTE modelling (TD-NLTE). 

Kneer (1980) showed that the relaxation timescale for 
the ionisation of hydrogen varies from 100 s to 1000 s in 
the middle to upper chromosphere. In the nineties, time- 
dependent hydrogen ionisation was first implemented in one- 
dimensional (ID) hydro dynamics simulations of the solar chro 
mosphere and corona (Carlsson & Stein 1992; Hansteen ll993 



ICarlsson & Stem! l2002h . ICarlsson & Stem! d2002l hereafter 
CS2002) gave a detailed analysis of dynamic hydrogen ionisa- 
tion in ID, where they included a t ransition region and a corona 
in their model. More recently, iRammacher & Ulmschneider 



(2003) also implemented a time-dependent treatment of mag- 
nesium in a ID simulation. 

The next step is obviously the inclusion of time-dependent 
hydrogen ionisation into 2D and 3D simulations. For reasons 
of numerical stability and current computer limitations, a full 
time-dependent NLTE treatment of hydrogen in the solar at- 
mosphere is not yet possible in two- and three-dimensional 
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models. In the ID case a simplified treatment using fixed ra- 
diative rates g ives a good approximation to the full solution 
( Sollum| [T999). In this paper, we report on the implementa- 



tion of this metho d in the 2D/3D radia tion hydrodynamics code 
C0 5 BOLD code dFrevtag et alJl2002h . It yields electron densi- 
ties that are more realistic than in the original simulation, as- 
suming LTE. The electron densities can now directly be used 
for a more realistic synthesis of chromospheric spectral lines 
and continua, producing much more reliable results compared 
to earlier calculations based on instantaneous statistical equi- 
librium. 

We give details on hydrogen ionisation in Sect. [2] and the 
numerical method in Sect. [3] The results of our numerical sim- 
ulation are described in Sect. |4] and finally discussed in Sect. [5] 

2. Hydrogen ionisation in the solar chromosphere 

In this section we recapitulate the details of chromospheric hy- 
drogen ionisation as found by CS2002. 

Statistical equilibrium does not hold under chromospheric 
conditions because the dynamical time scale of the chromo- 
sphere is much shorter than the relaxation timescale of hydro- 
gen level populations. The timescale in cool inter- shock regions 
was found to be of the order of 10 4 s. Within shocks the re- 
laxation timescale is typically 50 s, owing to the higher tem- 
peratures and densities. This is still sufficiently long to prevent 
equilibration of the hydrogen populations within a shock. Thus, 
the time-dependent rate equations need to be solved to obtain 
the correct hydrogen level populations. Up to just below the 
transition region Lyman a may be regarded as being in detailed 
balance for all practical purposes. Lyman continuum ionisation 
is negligible. To quote CS2002 literally: "the dominant hydro- 
gen ionization process is photoionization from the second level 
[...]. There is a very rapid equilibration of the second and all 
higher levels with the continuum, with a slow collisional leak- 
age of electrons from the ground state to the second level or 
visa versa, depending on whether hydrogen is ionizing or re- 
combining." Because the ionisation in shocks is much faster 
than the post-shock recombination, the ionisation degree and 
electron density tend to represent their shock value, also in the 
low temperature inter- shock regions. 

3. Method 

3.1. The approximate hydrogen radiative transport 

The full time-dependent NLTE radiation transport problem is 
as yet unsolvable in a 3D (M-)HD code. The demands on 
memory and processor time are very high, and the currently 
employed solution methods are not stable enough to handle 
strong shock waves in a 3D dynamical code. Simplifications 
are needed in order to compute the non-equilibrium hydrogen 
ionisation. We therefore adopt the concept of fixed radiative 
rates in order to avoid computationally expensive lambda iter- 
ation for evaluation of mean intensities. The new level popula- 
tions can be computed from the local gas density and tempera- 
ture, the imposed radiation field, and the level populations and 
electron densities from the previous time step. We so follow 



the example of Solium (1999) who showed that a simplified 
treatment of the radiative transport in hydrogen transitions, us- 
ing fixed radiative rates, gives reasonably accurate results when 
compared with both static and dynamical calculations in ID. 
He compared the full and simplified tr eatment both in t he static 
case of the VAL3C atmosphere ( Vernazza et al. 1981), and in 
dynamical sim ulations with the ID radiation hydrodynamics 
code RADYN dCarlsson & Steinl Il992h . He found that up to 
logTsoo = -5.6 the differences in the hydrogen level popu- 
lations between the full and the simplified treatment were at 
most a few percent. Above this height (just below the transition 
region in his model) large deviations occur because of the ne- 
glect of Lyman transitions in the simplified treatment. Below 
we will briefly outline his method and our extensions for 2D 
and 3D simulations. See lSolluml dl999|) for more details. 



3.1.1. Radiation field 

We assume that the radiation field in each transition, both 
bound-bound and bound-free, can be described by a formal ra- 
diation temperature r ra( j so that the angle-averaged intensity J v 
at frequency v is given by 
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where h, c, and k are Planck's constant, the speed of light, and 
Boltzmann's constant. The radiation temperature r ra( j is differ- 
ent for each transition but independent of frequency within a 
particular transition. However, r ra( j is allowed to vary over posi- 
tion and time. Its value is determined as follows: for each tran- 
sition the radiation temperature at the upper boundary of the 
model atmosphere r*°J is given as input, as derived by Sollum's 
detailed ID study. Then for each time step and each horizontal 
position we compute z C nt{x, y), the smallest height for which 



J V0 (T Q ) = 2/ V0 (O. 



(2) 



(3) 



Here T Q is the local gas temperature and vo the line centre fre- 
quency in the case of bound-bound transition and the ionisa- 
tion threshold frequency in the case of bound-free transitions. 
Below this height we set r ra( j = T e , i.e., we set the radia- 
tion temperature equal to the local gas temperature. Above this 
height, we set 

J V {Z) = B V {TZ) + tWefent)) " B V (T%)] ( -T^f 

Here m c (z) is the column mass at height z and H is a free pa- 
rameter that was determined by Solium. This means we have a 
decay of J v as a function of column mass. From this J v (z) we 
calculate T md (z) for z > z C rit- 

This procedure ensures the radiation temperature connects 
smoothly to the gas temperature in each column at z C rit- There 
may be strong horizontal gradients in the radiation temperature 
around this height if the gas temperature shows strong horizon- 
tal gradients as well. The most important radiative transition for 
hydrogen ionisation in the chromosphere is the B aimer contin- 
uum. In this transition the radiation temperature ranges from 
6250 K at z C rit to 5237 K at the top of the computational do- 
main. 
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3.1.2. Radiative bound-bound rate coefficients 



3.1.4. Collisional rate coefficients 



We restate the derivation of the radiative rate coefficients by 
Solluml dl999h . 



For lines, the assumption of a narrow line (J v /v is constant 
over the line profile) gives the following expression for the ra- 
diative excitation rate coefficient: 
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hv§m t c Jlu c 2 Q hv o/ kT md - l ' 
Here <x/ M (v) is the absorption cross section at frequency v, B[ u is 
the Einstein coefficient for radiative excitation, J VQ is the angle- 
averaged radiation field at the line centre frequency vo, and fi u 
is the oscillator strength. In addition we used the relation 



criu(y) dv = —B tu = fi u . 

Jo 47r m e c 

The bound-bound radiative deexcitation rate coefficient is 
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where gi and g u are the statistical weights of the lower and up- 
per level and A u \ and B u i are the Einstein coefficients for spon- 
taneous and stimulated deexcitation. 

3.1.3. Radiative bound-free rate coefficients 

The hydrogen bound-free excitations have a Kramers absorp- 
tion cross section 

(10) 
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where ao is the absorption cross-section at the ionisation edge 
frequency vo = (E c - Ei)/h, where E c and E[ are the energy 
of the continuum and the bound level. In this case the radiative 
excitation rate coefficient can be written as 
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where E\ is the first exponential integral and J y = B v (T r ^). The 
bound-free radiative deexcitation rate coefficient is 
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Using the Kramers cross-section and the radiation temperature 
this can be written as 
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Here n x , n c , and T e are the population of the lower level, the 
population of the continuum, and the local kinetic temperature, 
respectively. 



We take the collisional rate coefficients Cu from Ijohnsonl 



(I1972|) . who provides approximate cross sections for collisions 
of hydrogen with electrons. 

3.1.5. Electron densities 

As hydrogen is the most abundant element in the solar gas, 
ionisation of hydrogen severely affects the electron densities, 
which in turn enter into the collisional rate coefficients. LTE 
electron densities are a poor choice at chromospheric heights. 
Hence, we adopt a hybrid between LTE and TD-NLTE elec- 
tron densities. The solution of the rate equations in our TD- 
NLTE approach provides the population density of the contin- 
uum level H ii that is at the same time the contribution of hy- 
drogen to the electron number density. To this we add the LTE 
electron densities from the other elements by solving the Saha 
equation keeping the hydrogen ionization degree constant. For 
computational speed, we use a pre-calculated table that gives 
the electron density as a function of hydrogen ionisation de- 
gree, gas temperat ure, and mass density. Parts of the RH code 
( Uitenbroekl I200 1|) and partition functions by KurucJ3 were 
used for the creation of the table. 



3.1.6. Hydrogen level populations 

The equations and methods described above are used for solv- 
ing the level population evolution equations for hydrogen, 

+ V • ( B ,-v) = Yj n j p H ~ n i Tj P <i> (16) 

where ni, ni, and v are the total number of levels, the number 
density of level i, and the flow velocity. The second term on the 
left hand side of Eq. (IT6b represents changes in the population 
number due to advection by the hydrodynamic flow. The right 
hand side terms are rates into and out of level i by collisional 
and radiative transitions between the energy levels with Pj t and 
Pij being the rate coefficients for transitions from level j to i 
and vice versa. The rate coefficient is the sum of collisional 
and radiative rate coefficients: 



(17) 



(13) 3.2. Numerical implementation 



Numerically, Eq. (IT6l) is solved in two steps. In the first step, 
the populations from the previous time step are advected with 
the flow in the hydrodynamics computation. In the second step, 
the advected populations, mass density, and temperatures are 
used to calculate the new rate coefficients. Then the rate equa- 
tions (without the advection term) form a system of linear ordi- 
nary differential equations of first order, which is solved on the 
time interval between the previous and cu rrent time step using 
the DVODE package dBrown et al.lll989|) . The solution yields 
the level populations for the current time step. This system is 
solved with a relative accuracy of 10" 4 for each timestep. To en- 
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Fig. 1. Time evolution of the ionisation degree F in a column 
of the simulation for the first 30 min of the simulation. Upper 
panel: LTE. Lower panel: TD-NLTE. The TD-NLTE simu- 
lation was started with LTE populations. The chromosphere 
shows large temporal variation in ionisation degree in the LTE 
case. The ionisation degree in the TD-NLTE case is relatively 
constant after the first 5 min, when the first few shocks have 
passed. 



sure exact particle conservation we then scale the updated pop- 
ulation densities with the total hydrogen number density, which 
is directly calculated from the mass density and the mass frac- 
tion of hydrogen in the gas. This prevents cumulative buildup 
of small errors. The code outputs the electron densities and hy- 
drogen level populations for the LTE case and the TD-NLTE 
case for each grid cell. 

3.2.1. The hydrodynamic simulation 

We performed a three-dimensional numerical s imulation us- 
ing an upgraded version of the C O 5 BOLD code ( Frevta sTet alJ 
2002). The relaxed end model by Wedemeyer et al. ( 2004) was 
adopted as initial model. The computational grid consists of 
140 cells in each horizontal direction and 200 cells in the verti- 
cal direction. Horizontal resolution is 40 km, vertical resolution 
ranges from 46 km at the bottom to 12 km at the top of the com- 
putational domain. The horizontal extent is 5600 km. The verti- 
cal extent ranges from -1300 km below to 1783 km above aver- 
age Rosseland optical depth unity. Radiative transfer is treated 
in the grey approx imation. For more deta ils on the hydro- 
dynam ics code see IWedemeyerl d2003l) and I Wedemeyer et al.l 
( 2004). We let the simulation advance for 90 min of solar time. 
The first 30 min are reserved to ensure relaxation of the initial 



conditions, whereas the remaining 60 min are used for detailed 
analysis. 

3.2.2. The model atom 

We used a five-level-plus-continuum hydrogen model atom 
with collisional transitions from each level to each other level. 
Radiative transitions from the ground level are not allowed, i.e., 
all Lyman transitions are put in detailed balance. Lyman cool- 
ing or heating is thus not possible, but as our model does not 
include a transition zone no significant effect even in the upper 
layers of our model is to be expected. The bound-bound oscil- 
lator strengths come from I Johnson! (11972b whereas the bound- 
free radiative cross sections are from lSeato n (1960). 

4. Results 

Figure [T] shows the time evolution of the ionisation degree F 
along a column in the atmosphere. At t = min, the simulation 
is started with LTE values for the hydrogen level populations. 
The LTE panel (upper) shows the chromospheric shocks as red 
streaks of high ionisation in the upper part of the atmosphere, 
with cool, neutral gas in between. The typical time between 
successive shocks is of the order of 2-3 min. In contrast, the 
ionisation degree in the TD-NLTE case (lower panel) reaches 
a dynamical equilibrium state of fairly constant ionisation al- 
ready after 5 min, when the first shocks have passed. 

Figure [2] shows horizontal slices through the simulation af- 
ter 60 min. At z = 0.5 Mm the temperature and density show 
little structure, as it is in between the reversed g ranulation 
layer below (see lLeenaarts & Wedemeyer-Bohml2005b and the 
onset of strong shock formation in higher layers. The time- 
dependent ionisation degree is almost in LTE, showing only 
a somewhat higher ionisation degree in the cooler areas, where 
over-ionisation in the B aimer continuum occurs. At z = 1 Mm 
the situation is different. Owing to the steady decline in aver- 
age density, the waves that are excited by the convective mo- 
tions in the photosphere have steepened into shocks. Both the 
temperature and density fluctuations have increased. This has a 
profound effect on the LTE ionisation degree. Because the ion- 
isation degree is inversely proportional to the electron density, 
the hydrogen ionisation at equal temperature is much stronger 
than in deeper layers. In this particular snapshot the peak ion- 
isation degree at z = 1 Mm is 0.89 at T = 6813 K and the mini- 
mum is 3.5 x 10" 26 at T = 1871 K. In the time-dependent case 
the ionisation degree has less extreme values, the maximum is 
0.019 at 6443 K and the minimum is 7.8 x 10" 5 at 2495 K. Note 
also that the positions of the extrema of the TD-NLTE ionisa- 
tion degree do not at all coincide with the positions of the LTE 
extrema. The LTE ionisation degree depends on the local mass 
density and temperature only, whereas in the time-dependent 
case it also depends on the previous history of the atmosphere. 

At z = 1.5 Mm this effect becomes even more striking. The 
average density is lower, and temperature and density fluctua- 
tions are larger, resulting in even larger LTE ionisation degree 
fluctuations. In TD-NLTE the transition rates have become so 
small in the cool, non-shocked areas that there is almost no 
recombination between the passages of two shock waves. The 
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Fig. 2. Horizontal slices through a simulation snapshot at t = 60 min. Columns from left to right: gas temperature, logarithm of 
the mass density fluctuations log(p/(p)), LTE ionisation degree and TD-NLTE ionisation degree F. Top row: z = 0.5 Mm; middle 
row: z = 1 Mm; bottom row: z = 1.5 Mm. The LTE and TD-NLTE ionisation degree panels have the same color scale, which has 
been clipped at log F = -10. The black line indicates the position of the vertical cuts of Figs. [3] and ffl 



transition rates are higher in shock fronts as a direct result of 
higher temperature and density. Although this leads to an in- 
creased ionisation degree in shock fronts, it remains small com- 
pared to the LTE case. As a result the time-dependent ionisation 
degree at this height is fairly constant, in this snapshot varying 
only between 0.0014 and 0.021 with an average of about 0.008. 

Figure [3] shows some physical quantities along a vertical 
slice of the same snapshot as in Fig. [2] Panel a shows the 
temperature, with granulation near Mm, reversed granulation 
around 0.2 -0.3 Mm, and shocks above ~ 0.7 Mm. The mass 
density (panel b) decreases towards larger height, while its hor- 
izontal fluctuations increase. Panels c and d show the hydrogen 
ionisation degree F in LTE and TD-NLTE, respectively. The 
LTE case follows the increase of horizontal inhomogeneities 
with height, with high ionisation in the high-temperature shock 



waves. In contrast, the TD-NLTE case shows only small hori- 
zontal variation and the presence of shocks can hardly be dis- 
cerned. Finally, panels e and f show the electron density in LTE 
and TD-NLTE. The electron density is mainly set by the hydro- 
gen ionisation and - in first approximation - can be thought 
of as a multiplication of the density with F. Contributions 
from other elements are dominant only where F is lower than 
1 x 10" 4 , i.e., around 0.4 Mm and - in the LTE case - in the cool 
areas in the chromosphere. Again, in the time-dependent case 
the extrema and horizontal variations are smaller in TD-NLTE 
than in LTE. 

Figure |4] shows the departure coefficients for each level of 
our model atom in the same slice as Fig. [3] The departure coef- 
ficient is defined as the ratio of the population densities of a par- 
ticular atomic energy level in NLTE and in LTE, b x - rii/n^™. 
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Fig. 3. Vertical cut through the snapshot along the line indi- Fig. 4. Departure coefficients of our model hydrogen atom 

cated in Fig. [2] a: gas temperature; b: mass density; LTE (c) in the same cut as Fig. [3] for the continuum and level n=5 

and TD-NLTE (d) ionisation degree; LTE (e) and TD-NLTE (f) down to level n=l from top to bottom. The solution is nearly 

electron density. Note that the deepest layers of the model are in LTE from the bottom of the computational domain up to 

not shown. 0.3 Mm above the average Rosseland optical depth unity. The 

largest deviations occur in cool chromospheric regions in be- 
tween shocks, where, owing to the low temperature, the Saha- 
Below z = 0.3 Mm the departure coefficients for all levels Boltzmann equation predicts a very low occupation fraction for 
are close to 1, i.e., hydrogen is in or very near LTE. Above all excited levels, 
these heights deviations occur. The ground level (n = 1 , lowest 
panel) is very close to LTE for most of the chromosphere, ex- 
cept in strong shocks. In the shocks the ionisation lags behind slightly longer than the typical shock crossing time. The den- 
compared to LTE, resulting in a slight overpopulation of the sity of the ionised state H n (the continuum level) is correspond- 
ground level (which is the reservoir where ionised hydrogen ingly lower than in LTE. Due to the strong coupling to the H n 
atoms are coming from). The other levels are strongly colli- density,/? < 1 for all excited levels of neutral hydrogen in = 2 to 
sionally coupled with each other and show similar behaviour n = 5), too. In the cool intershock areas all levels except n = 1 are 
in the chromosphere. They are slightly underpopulated in the hugely overpopulated. This is because the Saha and Boltzmann 
shock area because the ionisation time scale is comparable or equations predict extremely low population fractions for ex- 
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Fig. 5. Histogram of the logarithmic ionisation degree as a 
function of height for the LTE (upper panel) and TD-NLTE 
(lower panel) case. The averages in LTE and TD-NLTE are 
plotted as red dashed and blue solid lines, respectively. All 
columns have been individually scaled to maximum contrast 
to enhance visibility. 



cited states and the continuum for temperatures and electron 
densities typical for the intershock regions. Thus, these enor- 
mous overpopulations simply illustrate the complete failure of 
Saha-Boltzmann equilibrium partitioning. 

Figure \5\ shows a height-dependent histogram of the ioni- 
sation degree. In LTE there is a strong bifurcation of the ion- 
isation degree as a result of the bimodal temperatu re distribu- 



tion i n our model chromosphere (cf. Fig. 7 of lWedemeyer et al 



2004)) . In TD-NLTE this bifurcation is not present. For both 
cases there is a minimum in average ionisation degree at z = 
0.3 Mm, roughly corresponding to the classical temperature 
minimum in ID static quiet- Sun models. Quite surprisingly, 
the average ionisation degree in the chromosphere is lower in 
the TD-NLTE case than in LTE. This is because the average 
is determined mainly by the high LTE ionisation degree in the 
shocks whereas in the TD-NLTE case the ionisation degree is 
smaller there. 



5. Discussion and conclusions 

Figure [6] shows a comparison of the average ionisation de- 
gree for our m odel, the detailed ID study of CS2002, and the 
FAL model C (iFontenla et al.lll993|) . The latter is a static and 
one-dimensional semi-empirical model, constructed under the 
assumption of statistical equilibrium, while the RADYN and 
CO 5 BOLD models are dynamic but differ in the number of spa- 
tial dimensions. All models nevertheless show the same quali- 
tative behaviour, with a minimum in ionisation degree between 
200 and 500 km, and a rise of the average ionisation degree 
with height in the chromosphere. The cause of this behaviour, 
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Fig. 6. Average TD-NLTE ionisation degree for our CO 5 BOLD 
model (solid), CS2002's RADYN model (dashed) and the sta- 
tistical equilibrium FAL model C (dotted). The zero point of 
the height scale is the average T500 = 1 height. 




Fig. 7. Vertical cut through the snapshot along the line indi- 
cated in Fig. [2 showing the relative contribution of other el- 
ements than hydrogen to the electron density. In the chromo- 
sphere the other elements contribute mostly in the high tem- 
perature shocks. 



however, is different. In the static FAL-C model, the ionisa- 
tion degree reflects a corresponding chromospheric tempera- 
ture rise. In contrast the rise in ionisation degree in the dynam- 
ical RADYN and CO 5 BOLD models is due to the temperature 
in shock fronts, which increa ses with height, witho ut need for 
a rise in average temperature dCarlsson & Stein|[l995 ). 

Qualitatively our CO 5 BOLD results are comparable to the 
detailed ID study of CS2002. We attribute the remaining dif- 
ferences in ionisation degree to the different dimensionality 
of the code (ID vs. 3D) and the resulting different temper- 
ature structure. ID simulations suffer fr om shock merging 



that l ead to too high shock temperatures dUlmschneider et al 
l2Q05h . Additionally, we do not include the effect of the time- 
dependent ionisation on the equation of stat e, causing too low 
shock temperatures in our CO 5 BOLD model ( Carlsson & Steinl 
Il992 ). Third, CS2002 include a transition region and corona in 
their model which influence the ionisation degree in the up- 
per chromosphere, both by Lyman radiation ionising the upper 
chromosphere and the varying height of their transition zone. 
Their transition zone is sometimes considerably lower than 
1.5 Mm, influencing the average ionisation degree. Fourth, it 
cannot be ascertained that the fixed radiative rates that we em- 
ploy, which are accurate in ID, are also accurate in 3D. This 
might affect our results as well. 
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Contributions to the electron density from other elements 
than hydrogen are still treated under the assumption of LTE. 
This will mostly affect the upper photosphere where the hydro- 
gen ionisation degree is so low that metals such as iron, mag- 
nesium, and silicon are the main electron donors. Hydrogen is 
the dominant donor in all other regions. However, in the chro- 
mospheric shocks the other elements contribute around 10% of 
the electrons (see Fig. [7]). Thus we expect a small error in the 
chromospheric elect ron density due to the LTE treatmen t of the 
other elements fsee lRammacher & Ulmschneiderl2003l for TD- 
NLTE effects on magnesium ionisation). For future simulations 
that take the transition region into account, a TD-NLTE treat- 
ment of helium might be of importance for the electron density 
as well, as the higher temperatures will lead to significant he- 
lium ionisation. 

All in all we conclude that our 3D simulation with time- 
dependent hydrogen ionisation produces reasonably realistic - 
if probably somewhat too low - results for the ionisation degree 
and electron density, given the level of necessary simplification 
and resulting increase in computational speed. 

With our method we can supply snapshots of 3D (M- 
)HD solar atmosphere simulations for detailed radiative trans- 
fer calculations containing time-dependent electron densities. 
Up until now most MHD simulations could only provide 
LTE electron densities, which are far from realistic in the 
chromosphere. The TD-NLTE electron densities are crucial 
in forward modelling of chromospheric diagnostics like the 
Ca I I H&K and infrared lines and (sub-)m illimetre continua 
(see iLeenaarts & Wedemeyer-Borrml 120061 for an example of 
the effect on TD-NLTE electron densities on the latter), as col- 
lisions between electrons and atoms provide the coupling of 
the state of the gas to the radiation field via the source func- 
tion, and affect the opacities of for example H free-free and H" 
bound-free and free-free transitions. 

As it is now, the code computes the ionisation degree and 
electron density from the state parameters defined by the hy- 
drodynamics simulation. Future improvement of the method 
should take into account the back-coupling of the ionisation 
to the equation of state. The deviation from instantaneous ion- 
isation equilibrium makes it impossible to use pre-computed 
lookup tables for the equation of state. In addition, the lookup 
tables that give the bin-averaged source function and opacity 
(computed in LTE as function of temperature and pressure) 
are then no longer consistent with the state of the gas. New 
methods have to be developed to handle this complex prob- 
lem. In particular, implementing such physics in a Riemann- 
solver as used in C0 5 BOLD is non-trivial. Nevertheless this 
work is currently in progress. Significant influence on chro- 
mospheric wave propagation and temperature structure in the 
chromosphere can be expected (cf . ICarlsson & Steinlll992 ). 

Another necessary improvement is the relaxation of the as- 
sumption of the fixed radiative rates. This is mainly important 
for the B aimer continuum, which is the driver of the ionisation. 
A radiative rate that is proportional to the mean grey intensity, 
or, in the case of multi-group opacity methods, the mean in- 
tensity in the continuum is an obvious choice. Such a method 
would take the variation in the average radiation field due to 
the granulation pattern into account. It will give more accurate 



results in the case of strong horizontal inhomogeneities in gas 
temperature and density in the upper photosphere, as it is the 
case in the presence of magnetic fields. 

Simulations for other stellar types than the Sun are pos- 
sible but involve re-adjustment of the radiation temperatures. 
This requires a similar detailed ID analysis and calibration as 
Solium performed for the Sun for each spectral type separately. 
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